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Abstract 


The  Digital  Orrery  has  been  used  to  perform  an  integration  of  the 
motion  of  the  outer  planets  for  845  million  years.  This  integration 
indicates  that  the  long-term  motion  of  the  planet  Pluto  is  chaotic. 
Nearby  trajectories  diverge  exponentially  with  an  e-folding  time  of 
only  about  20  million  years. 
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Abstract 

The  Digital  Orrery  has  been  used  to  perform  an  in¬ 
tegration  of  the  motion  of  the  outer  planets  for  845  million 
years.  This  integration  indicates  that  the  long-term  motion 
of  the  planet  Pluto  is  chaotic.  Nearby  trajectories  diverge 
exponentially  with  an  e-folding  time  of  only  about  20  mil¬ 
lion  years. 

The  determination  of  the  stability  of  the  Solar  System  is  one  of  the  old¬ 
est  problems  in  dynamical  astronomy,  but  despite  considerable  attention  all 
attempts  to  prove  the  stability  of  the  system  have  failed.  Arnold  has  shown 
that  a  large  proportion  of  possible  solar  systems  are  stable  if  the  masses,  and 
orbital  eccentricities  and  inclinations  of  the  planets  are  sufficiently  small  [1]. 
The  actual  Solar  System,  however,  does  not  meet  the  stringent  requirements 
of  the  proof.  Certainly,  the  great  age  of  the  Solar  System  suggests  a  high 
level  of  stability,  but  the  nature  of  the  long-term  motion  remains  undeter¬ 
mined.  The  apparent  analytical  complexity  of  the  problem  of  determining 
the  stability  of  the  solar  system  has  led  us  to  investigate  the  stability  by 
means  of  numerical  models.  We  have  investigated  the  long-term  stability  of 
the  Solar  System  through  an  845  million  year  numerical  integration  of  the 
outer  planets  with  the  Digital  Orrery  [2],  a  special  purpose  computer  for 
studying  planetary  motion. 

Our  earlier  integrations  [3]  showed  surprisingly  long  periods  in  the  mo¬ 
tion  of  Pluto,  of  order  137  million  years,  and  longer  periods  or  a  possible 
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secular  drift  in  the  inclination  of  Pluto.  The  motion  of  Pluto  is  involved  in  a 
surprisingly  large  number  of  resonances.  The  libration  of  the  resonant  argu¬ 
ment  permits  the  orbit  of  Pluto  to  cross  the  orbit  of  Neptune  without  having 
close  encounters  f4l.  Williams  and  Benson  [5]  also  found  that  the  argument 
of  perihelion  was  oscillating  about  tt/2.  Thus  the  perihelion  of  Pluto  remains 
far  from  the  line  of  intersection  of  the  orbital  planes  of  Pluto  and  Neptune. 
Williams  and  Benson  noted  that  two  other  resonance  conditions  were  pos¬ 
sibly  satisfied  by  the  motion  of  Pluto.  Using  the  frequencies  determined  in 
our  200  million  year  integrations,  we  find  that  these  two  resonances  are  sat¬ 
isfied  to  the  precision  of  our  measurements.  In  addition,  the  137  million  year 
components  that  we  observed  in  the  elements  of  Pluto  are  associated  with 
yet  another  commensurability.  In  Hamiltonian  systems  resonances  are  asso¬ 
ciated  with  both  enhanced  stability  and  instability.  Without  the  libration 
of  the  resonant  argument  or  the  argument  of  perihelion  a  Neptune  crossing 
orbit  would  be  highly  unstable.  Nevertheless,  resonances  are  also  associated 
with  most  prominent  zones  of  instability.  This  is  particularly  true  if  several 
resonances  strongly  affect  the  motion. 

The  very  long  periods  and  the  abundance  of  resonances  in  the  motion 
of  Pluto  compelled  us  to  pursue  longer  integrations.  Our  new  integration 
indicates  that  the  long-term  motion  of  the  planet  Pluto,  and  by  implication 
the  motion  of  the  whole  Solar  System,  is  weakly  unstable.  Exponential  di¬ 
vergence  of  nearby  trajectories  indicates  that  the  motion  of  Pluto  is  chaotic 
with  the  largest  Lyapunov  exponent  estimated  to  be  10_73yr-1. 

Our  Numerical  Experiment 

For  many  years  the  longest  direct  integration  of  the  outer  planets  was  the 
one  million  year  integration  of  Cohen,  Hubbard,  and  Oesterwinter  [6].  Re¬ 
cently  several  new  integrations  of  the  outer  planets  have  been  performed  [7,3,8] 
The  longest  was  our  set  of  200  million  year  integrations.  Our  new  integra¬ 
tion  is  significantly  longer  and  more  accurate  than  all  previously  reported 
long-term  integrations. 

In  our  new  integration  of  the  motion  of  the  outer  planets  the  masses 
and  initial  conditions  were  the  same  as  those  used  in  our  200  million  year 
integrations  of  the  outer  planets.  The  planet  Pluto  is  taken  to  be  a  zero- 
mass  test  particle.  We  continue  to  neglect  the  effects  of  the  inner  planets, 
the  mass  lost  by  the  Sun  due  to  electromagnetic  radiation  and  solar  wind, 
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and  general  relativity.  The  most  serious  limitation  of  our  integration  is  our 
ignorance  of  the  true  masses  and  initial  conditions.  We  believe  that  our  model 
is  close  enough  to  the  actual  solar  sytem  that  its  study  sheds  light  on  the 
question  of  stability  of  the  solar  system.  To  draw  more  rigorous  conclusions 
the  sensitivity  of  our  conclusions  to  the  uncertainties  in  masses  and  initial 
conditions,  and  to  unmodeled  effects  must  be  determined. 

Our  earlier  integrations  were  limited  to  100  million  years  forward  and 
backward  in  time  because  of  the  accumulation  of  error,  which  was  most 
seriously  manifested  in  an  accumulated  longitude  error  of  Jupiter  of  order 
50  degrees.  In  our  new  integrations  we  continue  to  use  the  12th  order  Stormer 
predictor,  but  a  judicious  choice  of  stepsize  has  reduced  the  numerical  errors 
by  several  orders  of  magnitude.  In  all  of  our  integrations  the  error  in  energy 
of  the  system  varies  nearly  linearly  with  time.  In  the  regime  where  neither 
roundoff  nor  truncation  error  is  dominant  the  slope  of  the  energy  error  as  a 
function  of  time  depends  on  stepsize  in  a  complicated  way.  For  some  stepsizes 
the  energy  error  has  a  positive  slope;  for  others  the  slope  is  negative.  This 
suggests  that  there  might  be  special  stepsizes  for  which  there  is  no  linear 
growth  of  energy  error.  By  a  series  of  numerical  experiments  we  indeed 
found  that  there  are  values  of  the  stepsize  where  the  slope  of  the  linear 
trend  of  energy  vanishes.  The  special  stepsizes  become  better  defined  as  the 
integration  interval  of  the  experiments  is  increased. 

We  chose  our  stepsize  on  the  basis  of  a  dozen  3  million  year  integrations, 
and  numerous  shorter  integrations.  For  our  new  long  integration  we  chose 
the  stepsize  to  be  32.7  days.  This  seemingly  innocuous  change  from  a  stepsize 
near  40  days  dramatically  reduces  the  slope  of  the  energy  error,  by  roughly 
three  orders  of  magnitude.  If  the  numerical  integration  were  truncation-error 
dominated,  for  which  the  accumulated  error  is  proportional  to  the  hn,  where 
h  is  the  stepsize  and  n  is  the  order  of  the  integrator,  then  this  reduction 
of  stepsize  would  only  improve  the  accumulated  error  by  about  a  factor  of 
about  10. 

In  our  new  integration  the  relative  energy  error  (energy  minus  initial 
energy  divided  by  the  magnitude  of  the  initial  energy)  over  845  million  years 
is  —2.6  x  10“ 10;  the  growth  of  the  relative  energy  error  is  still  very  nearly  linear 
with  a  slope  of  —3.0  x  10-19  per  year.  By  comparison  the  rate  of  growth  of 
the  relative  energy  error  in  our  200  million  year  integrations  was  1.8  x  10-16 
per  year.  The  errors  in  other  integrations  of  the  outer  Solar  System  were 
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comparable  to  the  errors  in  our  200  million  year  integrations.  The  rate  of 
growth  of  energy  error  in  the  one  million  year  integration  of  Cohen,  Hubbard, 
and  Oesterwinter  was  2.4  x  10~16  per  year.  For  the  6  million  year  integration 
of  Kinoshita  and  Nakai  the  relative  energy  error  was  approximately  5  x  10~16 
per  year.  For  the  LONGSTOP  integration  the  growth  of  relative  energy  (as 
defined  in  this  paper)  was  -2.5  x  10~16  per  year.  Thus  the  rate  of  growih 
of  energy  error  in  the  integration  reported  in  this  paper  is  smaller  than  all 
previous  long-term  integrations  of  the  outer  planets  by  a  factor  of  about  600. 

We  verified  that  this  improvement  in  energy  conservation  was  reflected  in 
a  corresponding  improvement  in  position  and  velocity  errors  by  integrating 
the  outer  planets  forward  3  million  years  and  then  backward  to  recover  the 
initial  conditions,  over  a  range  of  stepsizes.  For  the  chosen  stepsize  of  32.7 
days  the  error  in  recovering  the  initial  positions  of  each  of  the  planets  is  of 
order  10-5AU  or  about  1500km.  Note  that  Jupiter  has  in  this  time  travelled 
2.5  x  1015km. 

The  error  in  the  longitude  of  Jupiter  can  be  estimated  if  we  assume  that 
the  energy  error  is  mainly  in  the  the  orbit  of  Jupiter.  The  relative  energy 
error  is  proportional  to  the  relative  error  in  orbital  frequency  so  the  error  in 
longitude  is  proportional  to  the  integral  of  the  relative  energy  error:  AA  as 
tnAE(t)/ E,  where  n  is  the  mean  motion  of  Jupiter  and  t  is  the  time  of 
integration.  Since  the  energy  error  grows  linearly  with  time  the  position 
error  grows  with  the  square  of  the  time.  We  find  that  the  accumulated  error 
in  the  longitude  of  Jupiter  after  100  million  years  is  only  about  4  minutes  of 
arc.  This  is  to  be  compared  with  the  50  degree  accumulated  error  estimated 
for  our  200  million  year  integrations.  The  error  in  the  longitude  of  Jupiter 
after  the  full  845  million  years  is  about  5  degrees. 

We  have  directly  measured  the  integration  error  in  the  determination  of 
the  position  of  Pluto  by  integrating  forward  and  backward  over  intervals  as 
long  as  3  million  years  to  determine  how  well  we  can  reproduce  the  initial 
conditions.  Over  such  short  intervals  the  roundtrip  error  in  the  position  of 
Pluto  grows  as  a  power  of  the  time  with  an  exponent  near  2.  The  error  in 
position  is  approximately  1.3  x  10_J9t2  AU  (where  t  is  in  years).  This  growth 
of  error  is  almost  entirely  in  the  integration  of  Pluto’s  orbit;  the  roundtrip 
error  is  roughly  the  same  when  we  integrate  the  whole  system  and  when  we 
integrate  Pluto  in  the  field  of  the  Sun  only.  It  is  interesting  to  note  that 
in  the  integrations  using  the  32.7  day  stepsize  the  position  errors  in  all  the 
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planets  are  comparable.  Extrapolation  of  the  roundtrip  error  for  Pluto  over 
the  full  845  million  year  integration  gives  an  error  in  longitude  of  less  than 
10  minutes  of  arc. 


Features  of  the  orbital  elements  of  Pluto 

Our  new  integrations  confirm  and  extend  our  previous  report  of  extremely 
long  periods  in  the  orbital  elements  of  Pluto. 

For  the  first  450  million  years  of  our  integration  we  recorded  the  state 
of  the  system  every  499983  days  (about  1369  years)  of  simulated  time.  For 
the  second  400  million  years  we  sampled  16  times  less  frequently.  For  each 
planet  the  positions  and  velocities  were  converted  to  orbital  elements  relative 
to  the  center  of  mass  of  those  bodies  with  smaller  semimajor  axes  [9].  The 
elements  were  used  to  form  the  variables  h,  k ,  p,  and  q ,  defined  as 

h  —  e  sin  w  k  =  e  cos  w 

p  =  sin(r'/2)  sinO  q  =  sin(z/2)  cos  ft 

where  i  is  the  inclination,  e  is  the  eccentricity,  fl  is  the  longitude  of  the 
ascending  node,  and  va  is  the  longitude  of  the  perihelion.  The  variables  p 
and  q  specify  the  orientation  of  the  orbital  plane,  and  the  variables  h  and  k 
specify  the  eccentricity  and  the  orientation  of  the  orbit  in  the  orbital  plane. 
h,  k,  p,  and  q  are  natural  variables  for  looking  at  the  long-term  behavior  of 
planetary  systems. 

Figure  1  shows  the  variation  of  h  with  time.  The  largest  component  is 
the  3.7  million  year  regression  of  the  longitude  of  perihelion.  The  27  million 
year  component  we  previously  reported  is  clearly  visible,  as  is  the  137  million 
year  component.  The  quantity  p  (not  shown)  has  similar  features. 

Figure  2  shows  the  inclination  of  Pluto  as  a  function  of  time.  Besides 
the  major  3.8  million  year  component  we  can  clearly  discern  the  34  million 
year  component  we  previously  reported.  Although  there  is  no  continuing 
secular  decline  in  the  inclination,  there  is  a  component  with  a  period  near  150 
million  years  and  evidence  for  a  component  with  a  period  of  approximately 
600  million  years.  The  existence  of  significant  orbital  variations  with  such 
long  periods  is  quite  surprising.  For  quasiperiodic  trajectories  we  expect 
to  find  frequencies  which  are  low  order  combinations  of  a  few  fundamental 
frequencies  (one  per  degree  of  freedom).  The  natural  timescale  for  the  long 
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Figure  1:  The  orbital  element  h  —  e  sin  cc7  for  Pluto  over  845  million  years. 
On  this  scale  the  dominant  period  (the  3.7  million  year  circulation  of  the 
longitude  of  perihelion)  is  barely  resolved.  The  most  obvious  component  has 
a  period  of  137  million  years.  The  sampling  interval  was  increased  in  the 
second  half  of  our  integration. 
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term  evolution  of  a  quasiperiodic  planetary  systems  is  set  by  the  periods 
of  the  circulation  of  the  nodes  and  perihelia.  The  presence  of  significant 
components  with  periods  much  longer  than  these  indicates  strong  low-order 
resonances.  An  abundance  of  low-order  resonances  can  give  rise  to  chaotic 
behavior.  Pluto  is  involved  in  at  least  five  resonances,  a  sixth  is  possibly 
indicated  by  the  600  million  year  component  that  we  can  see  in  figure  2.  Each 
time  the  motion  of  Pluto  has  been  investigated  over  longer  intervals  longer 
periods  have  been  found.  This  is  suggestive  that  the  motion  is  chaotic. 

Deterministic  Chaotic  Behavior 

In  most  conservative  dynamical  systems  Newton’s  equations  have  both 
regular  solutions  and  chaotic  solutions.  For  some  initial  conditions  the  mo¬ 
tion  is  quasiperiodic;  for  others  the  motion  is  chaotic.  Chaotic  behavior  is 
distinguished  from  quasiperiodic  behavior  by  the  way  in  which  nearby  trajec¬ 
tories  diverge.  [10,11]  Nearby  quasiperiodic  trajectories  diverge  linearly  with 
time  on  average,  whereas  nearby  chaotic  trajectories  diverge  exponentially 
with  time.  Quasiperiodic  motion  can  be  reduced  to  motion  on  a  multidimen- 
cional  torus;  the  frequency  spectrum  of  quasiperiodic  motion  has  as  many 
independent  frequencies  as  degrees  of  freedom.  The  frequency  spectrum  of 
chaotic  motion  is  more  complicated,  usually  appearing  to  have  a  broad-band 
component. 

The  Lyapunov  exponents  measure  the  average  rates  of  exponential  diver¬ 
gence  of  nearby  orbits.  The  Lyapunov  exponents  are  limits  for  large  time  of 
the  quantity  7  =  ln(<f/<f0)/(t  — 10),  where  d  is  the  distance  between  the  trajec¬ 
tory  and  an  infinitesimally  nearby  test  trajectory,  and  t  is  the  time.  For  any 
particular  trajectory  of  an  n-dimensional  system  there  can  be  n  distinct  Lya¬ 
punov  exponents,  depending  on  the  phase  space  direction  from  the  reference 
trajectory  to  the  nearby  trajectory.  In  Hamiltonian  systems  the  Lyapunov 
exponents  are  paired;  for  each  non-negative  exponent  there  is  a  non-positive 
exponent  with  equal  magnitude.  Thus  an  m  degree- of- freedom  Hamiltonian 
system  can  have  at  most  m  positive  exponents.  For  chaotic  trajectories  the 
largest  Lyapunov  exponent  is  positive;  for  quasiperiodic  trajectories  all  of 
the  Lyapunov  exponents  are  zero. 

Lyapunov  exponents  can  be  estimated  from  the  time  evolution  of  the 
phase  space  distance  between  a  reference  trajectory  and  nearby  test  trajec¬ 
tories  [10,12].  The  most  straightforward  approach  is  to  simply  follow  the 
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trajectories  of  a  small  cloud  of  particles  started  with  nearby  initial  condi¬ 
tions.  With  a  sufficiently  long  integration  we  can  determine  if  the  distances 
between  the  particles  in  the  cloud  diverge  exponentially  or  linearly.  If  the 
divergence  is  exponential  then  for  each  pair  of  particles  in  the  cloud  we  ob¬ 
tain  an  estimate  of  the  largest  Lyapunov  exponent.  With  this  method  the 
trajectories  eventually  diverge  so  much  that  they  no  longer  sample  the  same 
neighborhood  of  the  phase  space.  We  could  fix  this  by  periodically  rescaling 
the  cloud  to  be  near  the  reference  trajectory,  but  we  can  even  more  directly 
study  the  behavior  of  trajectories  in  the  neighborhood  of  a  reference  trajec¬ 
tory  by  integrating  the  variational  equations  along  with  the  reference  tra¬ 
jectory.  In  particular,  let  y'  =  f(y)  be  an  autonomous  system  of  first-order 
ordinary  differential  equations  and  y (t)  be  the  reference  trajectory.  We  de¬ 
fine  a  phase-space  variational  trajectory  y  +  6y  and  note  that  6y  satisfies 
a  linear  system  of  first-order  ordinary  differential  equations  with  coefficients 
that  depend  on  y (£),  6y'  =  J  •  6y  where  the  elements  of  the  Jacobian  matrix 
are  Jtj  =  dfi/dy 

Lyapunov  Exponent  of  Pluto 

We  estimated  the  largest  Lyapunov  exponent  of  Pluto  by  both  the  vari¬ 
ational  and  the  phase-space  distance  methods  over  the  second  half  of  our 
845  million  year  run.  Figure  3  shows  the  divergence  of  the  logarithm  of  the 
phase-space  distance  in  a  representative  2- particle  experiment  and  the  growth 
of  the  logarithm  of  the  variational  phase-space  distance.  We  measured  the 
phase-space  distance  by  the  ordinary  Euclidean  norm  in  the  6-dimensional 
space  with  position  and  velocity  coordinates.  We  measured  position  in  AU 
and  velocity  in  AU/day.  Since  the  magnitude  of  the  velocity  in  these  units  is 
small  compared  to  the  magnitude  of  the  position,  the  phase  space  distance  is 
effectively  equivalent  to  the  positional  distance,  and  we  refer  to  phase  space 
distances  in  terms  of  AU.  For  both  traces  in  this  plot  the  average  growth 
is  linear,  indicating  exponential  divergence  of  nearby  trajectories  with  an  e- 
folding  time  of  approximately  20  million  years.  The  shapes  of  these  graphs 
are  remarkably  similar  until  the  two-particle  divergence  grows  to  about  an 
astronomical  unit,  verifying  that  the  motion  in  the  neighborhood  of  Pluto  is 
properly  represented.  A  more  conservative  representation  of  this  data  is  to 
plot  the  logarithm  of  7  versus  the  logarithm  of  time,  as  shown  in  figure  4. 
The  leveling  off  of  this  graph  indicates  a  positive  Lyapunov  exponent. 
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Figure  3:  The  exponential  divergence  of  nearby  trajectories  is  indicated  by 
the  average  linear  growth  of  the  logarithms  of  the  distance  measures  as  a 
function  of  time.  In  the  upper  trace  we  see  the  growth  of  the  variational 
distance  around  a  reference  trajectory.  In  the  lower  trace  we  see  how  two 
Plutos  diverge  with  time.  The  distance  saturates  near  80AU  when  the  Plu- 
tos  are  on  opposite  sides  of  the  Sun.  The  variational  method  of  studying 
neighboring  trajectories  does  not  have  the  problem  of  saturation.  Note  that 
the  two  methods  are  in  excellent  agreement  until  the  two- trajectory  method 
has  nearly  saturated. 
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Figure  4:  The  conventional  representation  of  the  Lyapunov  exponent  calcu¬ 
lation,  the  logarithm  of  7  versus  the  logarithm  of  time.  Convergence  to  a 
positive  exponent  is  indicated  by  a  leveling  off;  for  regular  trajectories  this 
plot  approaches  a  line  with  slope  minus  one. 


To  study  the  details  of  the  divergence  of  nearby  trajectories  we  expand 
the  early  portion  of  the  divergence  graph.  Figure  5  shows  the  logarithms 
of  the  distances  between  several  pairs  of  particles  in  the  cloud  of  Plutos  in 
our  integration  versus  the  logarithm  of  time.  The  separation  starts  out  as 
a  power  law  with  an  exponent  near  3/2.  Only  after  some  time  does  the 
exponential  take  off.  The  power  law  is  dominated  by  the  exponential  only 
after  the  rate  of  growth  of  the  exponential  exceeds  the  rate  of  growth  of 
the  power  law.  This  suggests  that  the  portion  of  the  divergence  of  nearby 
trajectories  that  results  only  from  the  numerical  error  fits  a  3/2  power  law 
and  that  this  error  “seeds”  the  exponential  divergence  that  is  the  hallmark  of 
chaos.  We  tested  this  hypothesis  by  integrating  a  cloud  of  test  particles  with 
the  orbital  elements  of  Pluto  in  the  field  of  the  Sun  alone.  The  divergence  of 
these  Kepler  “Plutos”  grows  as  3.16  x  10~17f3/2  AU.  This  is  identical  to  the 
initial  divergence  of  the  Plutos  in  the  complete  dynamical  system,  showing 
that  two-body  numerical  error  completely  accounts  for  the  initial  divergence. 

Only  the  second  half  of  the  integration  was  used  in  the  computation  of 
the  Lyapunov  exponents,  because  the  measurement  in  the  first  half  of  ou. 
integration  was  contaminated  by  over-vigorous  application  of  the  renormal¬ 
ization  method,  and  gave  a  Lyapunov  exponent  about  a  factor  of  four  too 
large.  The  renormalization  interval  was  only  275,000  years,  which  was  far  too 
small.  The  renormalization  interval  must  be  long  enough  that  the  divergence 
of  neighboring  trajectories  is  dominated  by  the  exponential  divergence  due 
to  the  sensitive  dependence  on  initial  conditions  rather  than  the  power  law 
divergence  due  to  the  accumulation  of  numerical  errors.  In  our  experiment 
the  renormalization  interval  should  have  been  greater  than  30  million  years. 
It  is  important  to  emphasize  that  the  variational  method  of  measuring  the 
Lyapunov  exponent  has  none  of  these  problems. 

Discussion 

Usually  the  measurement  of  a  positive  Lyapunov  exponent  provides  a 
confirmation  of  what  is  already  visible  to  the  eye,  chaotic  trajectories  usually 
look  irregular.  In  this  case  the  plots  of  Pluto’s  orbital  elements  do  not  look 
particularly  irregular  (see  figures  1  and  2).  The  irregularity  of  the  motion 
does  manifest  itself  in  the  power  spectra.  For  a  quasiperiodic  trajectory 
the  power  spectrum  of  any  orbital  element  is  composed  of  integral  linear 
combinations  of  fundamental  frequencies,  where  the  number  of  fundamental 
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Figure  5:  Common  logarithm  of  the  distance  between  several  pairs  of  Plutos, 
in  AU,  versus  the  common  logarithm  of  the  time,  in  years.  The  initial  seg¬ 
ment  of  the  graph  closely  fits  a  3/2  power  law  (dashed  line).  The  solid  line 
is  an  exponential  chosen  to  fit  the  long-time  divergence  of  Plutos.  The  expo¬ 
nential  growth  takes  over  when  its  slope  exceeds  the  slope  of  the  power-law. 


frequencies  is  equal  to  the  number  of  degrees  of  freedom.  The  power  spectrum 
of  chaotic  trajectory  is  more  complicated,  usually  appearing  to  have  some 
broad-band  component. 

A  portion  of  the  power  spectrum  of  Neptune’s  h  is  shown  in  figure  6.  For 
comparison  the  same  portion  of  the  power  spectrum  of  Pluto’s  h  is  shown 
in  figure  7.  Hanning  windows  have  been  used  to  reduce  spectral  leakage; 
only  the  densely  sampled  part  of  the  run  was  used  in  the  computation  of 
the  Fourier  transform.  The  spectrum  of  Neptune  is  quite  complicated  but 
there  is  no  evidence  that  it  is  not  a  line  spectrum.  On  the  other  hand  the 
spectrum  of  Pluto  does  appear  to  have  a  broad-band  component.  Note  that 
both  of  these  spectra  are  computed  from  the  same  integration  run,  using  the 
same  numerical  methods.  They  are  subject  to  the  same  error  processes,  so 
the  differences  we  see  are  dynamical  in  origin.  The  broad-band  character  of 
the  Pluto  spectrum  is  consistent  with  the  chaotic  nature  of  the  motion  as 
indicated  by  our  measurement  of  a  positive  Lyapunov  exponent. 

The  lack  of  obvious  irregularity  in  the  orbital  elements  of  Pluto  indicates 
that  the  portion  of  the  chaotic  zone  in  which  Pluto  is  currently  moving  is 
rather  small.  Since  the  global  structure  of  the  chaotic  zone  is  not  known  it  is 
not  possible  for  us  to  predict  whether  more  irregular  motions  are  likely.  If  the 
small  chaotic  zone  in  which  Pluto  is  found  connects  to  a  larger  chaotic  region 
relatively  sudden  transitions  can  be  made  to  more  irregular  motion.  This 
actually  occurs  for  the  motion  of  asteroids  near  the  3/1  Kirkwood  gap  [13]. 

On  the  other  hand,  the  fact  that  the  timescale  for  divergence  is  only  an 
order  of  magnitude  larger  than  the  fundamental  timescales  of  the  system 
indicates  that  the  chaotic  behavior  is  robust.  It  is  not  a  narrow  chaotic  zone 
associated  with  a  high  order  resonance.  Even  though  we  do  not  know  the 
sensitivity  of  the  observed  chaotic  behavior  to  the  uncertainties  in  parameters 
and  initial  conditions,  and  unmodelled  effects  the  large  Lyapunov  exponent 
suggests  that  the  chaotic  behavior  is  characteristic  of  a  range  of  solar  systems 
including  the  actual  solar  system. 

Conclusion 

Our  numerical  model  indicates  that  the  motion  of  Pluto  is  chaotic.  The 
largest  Lyapunov  exponent  is  about  10-7,3yr-1.  Thus  the  e-folding  time  for 
the  divergence  of  trajectories  is  about  20  million  years.  This  is  a  remarkably 
short  e-folding  time,  considering  the  age  of  the  Solar  System. 


Frequency  [10  cycles/day] 


Figure  6:  A  portion  of  the  power  spectrum  of  Neptune’s  h.  The  spectrum  is 
quite  complicated  but  there  is  no  evidence  of  a  broadband  component;  the 
spectrum  is  consistent  with  a  line  spectrum. 
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